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SYMBOLS 


a sound speed (includes sound speed perturbation due to wave passage) 
a Q ambient sound speed 

A ray tube area as cut by the wavefront 

Afo ray tube area as cut by a horizontal plane 

c n speed that a wave propagates normal to itself, given by a Q + V Q • N 
c Q defined as c n ! cos 6 

Cj atmospheric quantity defined by equation (34) 

C 2 atmospheric quantity defined by equation (35) 

F{%) an F function 

A Fj discontinuity of F(|) corresponding to A 

m i slope of waveform segment i 

m® initial value of 

M Mach number 

N wavefront unit normal 

p pressure perturbation due to wave passage 

p Q ambient pressure 

Apj pressure rise across shock at the juncture of waveform segments i and i — 1 
A pP initial value of A Pj 

r for a conical wave, the distance from the flight path of the projectile that is generating the 
wave; for a spherical wave, the distance from the point where the wave originated 

r Q initial value of r 

S for a plane wave, the distance of propagation 
t time 

iii 



At propagation time required for the wave to traverse the distance between points on the ray 
path 

T time associated with a waveform point; also, a time parameter defined by equation (39) 

T e defined by equation (25) 

T f defined by equation (26) 

c 2 

u fluid particle speed due to wave passage (positive in the direction of wave propagation) 
V Q wind velocity 

z distance above ground 

Zj initial altitude, where the pressure signature is known 
z 2 altitude where the pressure signature is desired 
Cj defined by sketch (c) and equation (20) 

e 2 defined by sketch (c) and equation (20) (if in eq. (20) is replaced by e 2 and A Fj is 
replaced by A Fj + j ) 

7 ratio of specific heats 

time duration of waveform segment i 

\° initial value of X ; - 

p Q ambient density 

t age variable 

9 angle between N and horizontal plane 

£ phase, a variable used to identify or “tag” points on a waveform 
Subscript Notation for Derivatives 

( ) z = 3( )/9z 

( )$ = 3 ( )/ 9 € 

( )^ z = a H )/9za^ 


IV 


EXTRAPOLATION OF SONIC BOOM PRESSURE SIGNATURES 


BY THE WAVEFORM PARAMETER METHOD 
Charles L. Thomas 


Ames Research Center 


SUMMARY 


The waveform parameter method of sonic boom extrapolation is derived and shown to be 
equivalent to the F-function method. A computer program based upon the waveform parameter 
method is presented and discussed, with a sample case demonstrating program input and output. 


INTRODUCTION 


Sonic boom pressure waves created by high-flying supersonic aircraft propagate over ex- 
tremely large distances through an atmosphere with temperature, pressure, and wind gradients 
before reaching the ground. The atmospheric gradients have a strong effect on both the wave am- 
plitude and the nonlinear waveform distortion. Two sonic boom extrapolation methods that ac- 
count for these atmospheric effects are the F-function method (ref. 1) and the waveform parameter 
method (ref. 2). These methods are based on the same fundamental concepts; both rely on a result 
from geometric acoustics for the wave amplitude, and both utilize isentropic wave theory to ac- 
count for nonlinear waveform distortion. They differ, however, in the manner in which they ac- 
count for nonlinear effects. The F-function method uses an age variable, calculated from the 
atmospheric properties, to distort the pressure signature. Shock waves are then located from the 
distorted signature by means of an area-balancing criterion. No shock-finding criterion is required 
in the waveform parameter method. Instead, a set of waveform parameters is defined that com- 
pletely describes the waveform, and equations are obtained for the time rates of change of these 
parameters. These equations, which can be integrated over small time increments, are derived from 
the same basic equations used in the F-function method. No additional assumptions are required, 
and therefore the two methods are mathematically equivalent. The waveform parameter method 
appears to provide a more suitable approach for automatic computation because the necessity of 
using the area-balancing technique for locating shocks is eliminated. Also, in the special cases of 
plane, conical, and spherical pressure waves propagating in a uniform medium, closed-form expres- 
sions for the waveform parameters as functions of propagation distance can be obtained. 



ANALYSIS AND DISCUSSION OF RESULTS 


Both the F-function method and the waveform parameter method depend, for waveform 
amplitude, on a result derived from geometric acoustics— conservation of the Blokhintsev energy 
invariant along ray tubes (refs. 3 and 4). This conservation principle may be stated as 


p fnJi=Ftt) 0) 

V p o a o 

The quantities c n . A, p Q , and a Q are functions only of z (altitude) and therefore vary along the 
ray tube. The acoustic pressure p is a function of both z and phase £. If the sonic boom pres- 
sure wave were truly an acoustic wave, the above expression would be sufficient for calculation of 
the sonic boom waveform at any point along the ray tube. However, because of the extremely large 
propagation distances involved, and because near the aircraft the wave is not extremely weak, the 
actual nonlinearity of pressure wave propagation causes a distortion of the waveform that must be 
taken into account. Both the F-function method and the waveform parameter method account for 
this nonlinear distortion by assuming that the speed of propagation of any point on the waveform 
is given by isentropic wave theory. The propagation speed of a point on a pressure wave is equal to 
the value of u + a for the point. If the wave is assumed isentropic and of small amplitude, the 
propagation speed is given by 


u + a = 




( 2 ) 


or, in terms of the function F(£) 


u+a=a Q +ll J [j°Zm (3) 

2 \j p Q Ac n 

Let the pressure signature at some point along the ray tube be described by the pressure p 
versus time T plot that would be produced by a microphone located at that point. Also, let 
T - 0 correspond to a zero disturbance point on the waveform as shown by the example signature 
of sketch (a). Then the amount of nonlinear waveform distortion that occurs in a propagation time 
increment dt is given by 


dT-- a - o_zhl* dl 


(4) 


If the altitude z of the wave is measured from the T = 0 point on the waveform, then 


dz 

dt 


~a Q sin 0 


(5) 


2 


200 r 



Sketch (a) 


and, using equation (3), the total nonlinear distortion that occurs during propagation from the 
initial altitude z x to the altitude of interest z 2 is given by 


A7X£) = j 


z 2 dz 

l / p Q a 0 c n 4 A sin 2 d 


( 6 ) 


The phase £ is as yet undefined quantitatively. If £ for a point on the waveform is taken 
to be the value of T corresponding to the same point on the initial signature — that is, the signa- 
ture at Zj — then the required expression for 7T£) is 

T=i-rm) ( 7 ) 


where the “age variable” r is calculated from 


T = 


7+1 

~2~ 



2 dz 

'ff> 0 a o c n A sin 2 d 


( 8 ) 


It is pointed out by Hayes (ref. 1) that in a horizontally stratified, steady atmosphere 
c n ~ c o cos ® where c Q is constant along ray paths. It is possible to define the ray tube area as the 
area cut by a horizontal plane, rather than that cut by the wavefront. The relationship between the 
two areas can easily be shown to be 


A, =_£zl A — 
h a Q sin0 


(9) 


The expression for the age variable in terms of is then 


T = - H 1 


r + 1 r 

2 c n ™ J 

o * 7 


dz 


Z ^\Po a o A h si ° 3 6 cos3 °) 


Vx 


( 10 ) 


3 


In cases where sonic cutoff occurs, or nearly occurs, the use of A ^ instead of A presents some 
difficulty due to sin 6 going to zero. 

The above expression for r is identical to that obtained in reference 1 (eq. (46)), except for 
the exponent of c Q , which in reference 1 is 1 instead of 3/2. This difference in exponent is due to 
a factor of c 0 ~' /z in the definition of F(£) used in reference 1, where F(£) is denoted VgtZ). 
Since c Q is constant along ray tubes, including c Q in the definition of F(£) is permissible. 

The above equations are the fundamental equations of the F-function method. The F func- 
tion F(£) is first determined using a near-field signature and equation (1). The age variable cor- 
responding to altitude z 2 is then calculated from equation (8) or equation (10). Then, using 
equations (1) and (7), one can plot p versus T at altitude z 2 - The result will normally look as 
shown in sketch (b). It is seen that the signature (indicated by the solid line) is multivalued in p. 
The multivalued regions are eliminated by locating shock waves (indicated by dashed lines) accord- 
ing to the area-balancing rule. The area-balancing rule is a consequence of the fact that the propa- 
gation speed of a weak shock is equal to the average of the propagation speeds (u + a) just ahead 
of and just behind the shock. 



Sketch (b) 


In the waveform parameter approach, the waveform is approximated by an arbitrary number 
of linear segments. For each segment the waveform parameters m A p^, and X ? - are defined as fol- 
lows: m ; - is the slope bp/dT of segment i, which may be positive or negative; A p^ is the pres- 
sure rise across the shock at the juncture of segments i and /— 1. Often there will be no shock at 
the juncture, in which case A p^ is zero. Finally, X z - is the time duration AT of segment i. A com- 
pletely general waveform can be described using these waveform parameters. 

Consider first the waveform parameter mp which is defined by 



Then 


dm i- P Jth. di + P _hdz ( 12 ) 

dt T 2 dt ~ dt 

\ T it 

The first term on the right-hand side of equation (12) is zero for a linear, nonplane wave, whereas 
the second term is zero for a nonlinear, plane wave propagating in a uniform, stationary medium. 
This means that the rate of change of the waveform parameter m ^ can be obtained by super- * 
position of the rate of change assuming the wave propagates as a linear, nonplane wave and the 
rate of change assuming the wave propagates as a nonlinear, plane wave. This concept of super- 
position of waveform parameter rates of change was first introduced, without proof, in reference 2. 
Equation (12) can be written as 


Then, from 


dt dt p £ 9 1 


(13) 


p *t 


\Pn a 


O 0 Z7' 


c n A 


F/iZ) 


(14) 


7>z, = ~ F i (£) t z 


(15) 


7+1 

T z 


2 V p o a o c n A sir * 2 6 


(16) 


and equation (5), the time rate of change of the waveform parameter m j is given by 


1 %_ m 2 + Ul d ^^l_ d _£o_2, ( K_\dA\ m . 

dt 27 p o c n 1 2 \ o dt p o dt c n dt Adt I 1 


(17) 


where the fractional time rates of change of the ambient properties a Q , p Q , and c n and of the ray 
tube area A are understood to be the rates of change as seen by an observer moving down the ray 
tube with the wave. 
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The waveform parameter A p i can be expressed in terms of the jump in the F function 
across the shock as 


Then 


Ap/ = 



AF t 


(18) 


dApj _ Ap i dAFj 
dt A Fj dt 



+ A Pi ~ — , ~JT 


(19) 


The first term on the right-hand side of this equation is zero for a linear, nonplane wave, whereas 
the second term is zero for a nonlinear, plane wave propagating in a uniform, stationary medium. 
Therefore, the concept of superposition of waveform parameter rates of change is also valid for the 
parameter A p i The fractional time rate of change of A F i can be determined as follows. Sketch (c) 
shows a typical signature at some point along a ray tube. The dashed lines are “equalizing lines” 
representing the shock waves at some time increment dt later. Using equation (3), the relationship 
between Cj and dt is 


2e i c n 



A F t dt 


( 20 ) 



Sketch (c) 


Since from the area-balancing rule 


dAF t = e x 
dt dt 


Im i + m 




( 21 ) 
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we have 


_L d ^i 
A Ft dt 


7 + 1 a Q / 

4 7 P 0 c„\ 


+ m i 


-) 


( 22 ) 


and therefore the time rate of change of the waveform parameter A is 


+ m /-«) + 4 -4 #W,- 


* 2 wz 0 cfr p Q <2/ c w A dt J ' 


(23) 


The waveform parameter X ; - is defined as 


where 7^- and T ?+1 are the times corresponding to the shocks at each end of the waveform seg- 
ment, as shown in sketch (c). Let 


= Tj + e, 

(25) 

= T i + 1 ~ e 2 

(26) 


Since in time dt both e t and e 2 go to zero (by definition of Cj and e 2 ),then 


^}=d( T - T \- e -±^- e ' 
dt dt\ e 2 e ij dt 


(27) 


Let F c and F e be the values of the F function F(£) at T e and T p , respectively. Then, 
e i e 2 fc i e 2 

from equation (7) 




(28) 


= -™;V 


c n dr 


Pn a 


3 dt 


oo 


(29) 
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Thus, from equations (5), (16), (18), and (20) the time rate of change of the waveform parameter 
X z - is 


d\ 

~dt 



7+1 a o 
Po c n 


m i\ 


(30) 


We now have three ordinary, first-order, coupled differential equations for the waveform 
. parameters mp A P p and \p which completely describe the deformation of the waveform. These 
equations can be written more concisely as 


am : 0 

— = C, m , 2 + C, m f 
dt 1 1 21 


^Hr = i c i A Pi ( m i + m i~i) + C 2 A p i 

-\ c i ( A Pi + A Pi+ 1) - c i m i\ 

where 

^ - 7 + 1 a o 

1 ~W Pcfn 


r - 1 / 3 d _% + 1 d ±o _ 2 dc n _ 1 dA \ 
2 2\a 0 dt p Q dt c n dt A dt / 


(31) 

(32) 

(33) 

(34) 

(35) 


The above equations differ from corresponding equations of reference 2 because the waveform 
parameters used here refer to a pressure versus time signature, while those of reference 2 refer to a 
pressure versus distance signature (the distance being measured normal to the wavefront). 

In the general case of a wave with arbitrary wavefront shape propagating in a nonuniform at- 
mosphere with winds the quantities Cj and C 2 vary along the ray path. However, if these quanti- 
ties are assumed to be constant over small time increments, the above waveform deformation equa- 
tions can be integrated to obtain the following solutions. 


m t =- 


o C i At 
m^e 1 

1 -C^mfT 


(36) 


A Pi = 


[(' 


A Pi °e ~ 2 


C.At 


(l -C,<, 



(37) 


8 




Cj m 


i° T ) 




I l i - c i < i ? ~ _ A 

i y 1 -C^mpT J 


*Pi + 1 /l l - c i<iT \ 
<-<1 \Al i-qmfr y 


where 


T = 


e C 2 Ar - 1 


( 38 ) 


(39) 


The above solutions for A and X ; - are used to determine the values of the waveform param- 
eters at some point on the ray path, given m° , A p-°, and , which are the values of the param- 
eters at the preceding point. The time increment At is the propagation time required for the wave 
to traverse the distance between points on the ray path. To use the above solutions for the wave- 
form parameters one must first calculate the values of C- L and C~ at many points along the ray 
path of interest. The ray path can be calculated by the method outlined in reference 2. When the 
ray path is known, the fractional rates of change of the ambient properties a Q , p Q , and c n can be 
determined at many points along the ray path. The fractional rate of change of ray tube area is also 
needed at each point along the ray path. Ray tube area may be calculated by theoretical methods, 
such as that of reference 1 , or can be determined directly by independently calculating each of the 
four rays that bound the ray tube. Once the values of C l and C 2 are determined, the above solu- 
tions for the waveform parameters can be used to calculate the waveform at any point of interest 
along the ray path if the waveform near the aircraft is known. 

The above expression for X ; - (eq. (38)) cannot be used when m° = w/Lj or when 
m° = mf + i . In the case when mf* = or, in general, when 


c i( m i°- m ?-i) T 

1 -C x mfT 


< 0.001 


(40) 


we can use the approximation 


T _ 

W 1 -C x mp 



Apf CJ 

2 1 — Cj m^T 


(41) 


A similar expression can be obtained when mf = 
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As the wave propagates down the ray path, shocks often will coalesce and compression regions 
will steepen to form new shocks. When this occurs one or more of the X t - will go to zero. When one 
of the \j does go to zero somewhere between two points on the ray path, the expression for X z - 

given by equation (38) can be used to determine the value of At at which it goes to zero. The 
waveform parameters should then be incremented using this value of At to determine the wave- 
form at the point where the waveform segment goes to zero. The waveform parameters must then 
* be redefined so that the new waveform (which has one less segment and possibly one more shock) 
is represented correctly. 

Extrapolations of plane, conical, and spherical pressure waves propagating in a uniform 
medium can be done easily by hand if the waveform is of a very simple shape. Solutions of the 
waveform deformation equations for plane, conical, and spherical waves are presented in appendix A. 

A sonic boom extrapolation program, based on the waveform parameter method described 
above, is presented and discussed in appendix B. The program includes the effects of aircraft ac- 
celerations, aircraft flight path angle, and atmospheric temperature, pressure, and wind gradients. 
The program has been checked using the well-established Hayes program (ref. 1) and has been 
found to give nearly identical results for all types of flight conditions. 


CONCLUDING REMARKS 


The waveform parameter method for extrapolation of weak pressure waves has been shown to 
be valid, within the same limitations as the F-function theory. The waveform parameter method 
appears to provide a more suitable approach for automatic computation because the necessity of 
using the area-balancing technique for locating shocks is eliminated. The waveform parameter 
method also makes possible closed-form solutions for plane, conical, and spherical waves propa- 
gating in uniform media. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif. 94035, February 7, 1972 
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APPENDIX A 


SOLUTIONS OF THE WAVEFORM DEFORMATION EQUATIONS 
FOR PLANE, CONICAL, AND SPHERICAL WAVES 


Plane, conical, and spherical pressure waves propagating in a uniform, stationary gas can be 
treated very simply by the waveform parameter method. The waveform deformation equations 
(31), (32), and (33) can be solved in closed form in these special cases. The solutions are found to 
be the following. 





APj = A pP 




where the quantities Q 1 and Q 2 are given in the following table. 


Qi 


Q 2 


Plane wave 


5 


Conical wave 



Spherical wave 


r 
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For a plane wave, S is the distance of propagation. For a conical wave, r is the distance from 
the flight trajectory of the projectile that is generating the wave, and M is the Mach number of the 
projectile. For a spherical wave, r is the distance from the point where the wave originated. In 
both cases, r Q is the value of r corresponding to the initial pressure signature. 

If the value of Xj computed from the above equation is negative, the equation for \j must 
be used to determine the value of r or S at which becomes zero. The initial signature should 
then be redefined to be the signature at this value of r or S. 

It should be remembered that the above solutions for the waveform parameters assume a 
weak, isentropic wave. If p/p 0 < 0.1 over the entire signature, this assumption is justifiable. If the 
wave is conical or spherical and the initial signature is near the point where the wave originated, the 
waveform amplitude will attenuate rapidly due to a large fractional rate of change of ray tube area. 
In these cases, the maximum value of p/p 0 can be somewhat larger than 0.1 and the above expres- 
sions for the waveform parameters will still be applicable. 
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APPENDIX B 


A COMPUTER PROGRAM FOR SONIC BOOM EXTRAPOLATIONS 


The sonic boom extrapolation program presented here is based on the waveform parameter 
method described in this report. This computer program has essentially the same capability as the 
extrapolation program of reference 1. Furthermore, both programs have the same limitations, 
namely, those of geometric acoustics and isentropic wave theory. The advantage of the waveform 
parameter approach is that it provides a more appropriate base for automatic computation since no 
shock -finding technique is needed. As a result, the program presented here is considerably shorter 
and simpler than the program of reference 1 . Other advantages of the program presented here are 
the following. No F-function calculation is required, thus simplifying extrapolations of experimen- 
tal signatures. No problems arise in extrapolating near field signatures corresponding to multivalued 
F functions. The atmospheric pressure variation with altitude is computed within the program, and 
therefore does not have to be input. Aircraft accelerations are more conveniently specified. Sonic 
boom “footprints” are more easily determined since both aircraft location and ground-ray inter- 
sections are specified in terms of geographical longitude and latitude. The sonic cutoff angle is 
computed, thus eliminating timely trial-and-error determinations of this angle. Signatures at more 
than one altitude can be output on each run. Also, signatures can be output in terms of A p (psf), 
A p/p, x (ft), or t (msec). Finally, the effect of varying aircraft length can be easily determined, 
without requiring recalculation of the aircraft F function for each aircraft length. Because of these 
advantages the extrapolation program presented here is considerably more convenient to use than 
the program presented in reference 1 . 


In this extrapolation program ray path calculations are performed by the method outlined in 
reference 2. Ray tube areas are determined directly by independently calculating each of the four 
rays that bound the ray tube. 

Input to the program consists primarily of aircraft flight conditions, atmospheric properties, 
and near-field signature data. Flight conditions include Mach number, altitude, flight path angle, 
and aircraft accelerations. The three components of aircraft acceleration are expressed in terms of 
M, 7 , and \j/, which are the time rates of change of Mach number, flight path angle, and heading, 
respectively. Atmospheric properties are input in the form of temperature and wind velocity pro- 
files. The atmospheric pressure variation with altitude is computed in the program from the input 
temperature profile, using the perfect gas law and the hydrostatic equation. Winds in both the 
northerly and easterly directions are input, allowing wind shear. However, vertical winds arid atmos- 
pheric turbulence are not accounted for in this program. The near-field signature data consist of the 
signature itself and its corresponding location relative to the aircraft. The signature is input in the 
form of p/p 0 (more commonly written as Ap/p) versus x, where x is a spatial coordinate meas- 
ured parallel to the aircraft velocity vector. The signature that is input should be consistent with 
the flight conditions specified. The flight conditions that affect the near-field signature are Mach 
number and lift coefficient. Therefore, the lift coefficient should first be estimated from the air- 
craft weight, flight altitude, Mach number, and aircraft accelerations. The near-field signature cor- 
responding to the Mach number and required lift coefficient can then be determined either experi- 
mentally by a wind tunnel test or by theoretical means. Normally, if near-field signatures are 
determined in the wind tunnel, signatures are obtained at several lift coefficients, or angles of 
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attack, and in several azimuthal planes of the model. A near-field signature corresponding to a 
specific lift coefficient and location relative to the model can then be estimated by interpolation. 

Program output consists of a listing of the input data and of output signatures at the desired 
altitudes. The location (latitude and longitude) where the output signature would be detected is also 
output so that sonic boom “footprints” can be determined. In the event that sonic cutoff occurs, 
the sonic cutoff altitude is output. Also, the minimum </> (azimuthal angle) for which sonic cutoff 
occurs is output. In the event that a ray tube area goes to zero, due to large aircraft accelerations, 
the program will output the altitude at which the ray tube area goes to zero and the location of the 
ground-ray intersection. The computation will then be terminated. Because the basic theory is one 
of geometric acoustics, the theory breaks down in regions where ray tube areas go to zero. All that 
can be said regarding the wave amplitude is that a “superboom” is likely to occur if the ray tube 
area goes to zero near the ground. Although the program cannot be used to predict the amplitude 
of a superboom, it can be used to predict where superbooms are likely to occur. 
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Program Input Variables 


TITLE 

MACH 

FLTALT 

MDOT 

PSIDOT 

GAMDOT 

FPA 

LONGP 

LATP 

HEAD 

POG 

NTEMP 

NWINDE 

NWINDN 

ZT(I) 

TO(I) 

ZWE(I) 

VOE(I) 

ZWN(I) 

VON(I) 


title card. 80 columns available, 
flight Mach number 
flight altitude above ground in feet, 
time rate of change of MACH (/sec) 
time rate of change of HEAD (deg/sec) 
time rate of change of FPA (deg/ sec) 
aircraft flight path angle (deg) 

aircraft longitude (deg). West longitudes are taken to be positive. East longitudes are 
negative. 

aircraft latitude (deg). North latitudes are positive. South latitudes are negative, 
aircraft heading in degrees from true north. 

atmospheric pressure on the ground in psf. Input POG = 1 if output signatures in 
terms of Ap/p are desired. Input POG = 0 for uniform atmosphere extrapola- 
tions, or for any case where constant atmospheric pressure is desired. 

number of atmospheric temperature values input. NTEMP > 2 

number of east wind values input. NWINDE > 2 

number of north wind values input. NWINDN > 2 

altitudes (in thousands of feet) at which atmospheric temperatures are given. Pro- 
gram logic requires ZT(NTEMP) > FLTALT. 

atmospheric temperatures in degrees Fahrenheit. The program assumes a linear var- 
iation of temperature between points. 

altitudes (in thousands of feet) at which east wind values are given. Program logic 
requires ZWE(NWINDE) > FLTALT. 

wind speeds in east direction (blowing east) in ft/sec. The program assumes a linear 
variation of wind speed between points. 

altitudes (in thousands of feet) at which north wind values are given. Program logic 
requires ZWN(NWINDN) > FLTALT. 

wind speeds in north direction (blowing north) in ft/sec. The program assumes a 
linear variation of wind speed between points. 
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NALT number of altitudes at which output signatures are desired. 

ALT(I) the altitudes (feet) at which signatures are desired, in descending order. 

OCODE output code. OCODE = 1 for output signatures of pressure versus 

time (msec). OCODE = 2 for output signatures of pressure versus 
distance (feet), as measured along a line parallel to the horizontal component of the 
aircraft velocity (relative to the ground). 

REFL sonic boom reflection factor. Signatures on the ground are multiplied by REFL. 

ROVERL radial distance (nondimensionalized by the aircraft reference length) from the flight 

path corresponding to the input pressure signature. 

PHI the lateral angle in degrees corresponding to the input pressure signature and aircraft 

bank angle. PHI is measured from the vertical plane passing through the aircraft 
velocity vector. Positive values of PHI correspond to rays that start out to the left of 
vertical, as seen from behind the aircraft. 


Aircraft in left turn 
as seen from behind 

Plane through / I 

aircraft velocity vector^ / ] 


Near field signature 
at this point is input 


Sketch (d) 

NX number of values input for X(I) and DPP(I) 

X(I) The near field pressure signature is input as a Ap/p versus x signature, where x 

is a spatial coordinate measured parallel to the aircraft velocity vector. X(I) is a 
table of x values for the near field signature. X(I) can be input in any unit of 
length. 

DPP(I) the value of Ap/p corresponding to X(I). The program assumes a linear variation 

of Ap/p between points. 
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AL aircraft reference length in feet. 

ML model reference length in the same unit of length as used for X(I). 
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PROGRAM INPUT - SAMPLE CASE 


SPACE 

SHUTTLE 

BOOSTER 

REENTRY 

STANDARD 

ATMOSPHERE 

EASTERLY WINDS 

1.20 

50400. 









-.0197 

-.359 

1.013 

-12.75 







119.88 

2116.2 

27.60 

356.5 








7 

8 

2 








0 . 

36.2 

65.8 

105.5 

155.5 

172.0 

202.0 




8 9.0 

-69.7 

-69.7 

-48.1 

27.5 

27.5 

-4.8 




0 . 

30. 

40. 

45. 

56. 

65. 

83. 

110. 



5. 

68. 

84. 

79. 

36. 

19. 

16. 

34. 



0 . 

202. 









0 . 

0 . 









2 










30000. 

i 

0 . 









l 

1.9 

5.30 

47 . 
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6.8 

7.1 

7.4 

7.6 

7.8 

8.7 

9.1 

9.4 

9.6 

10.0 

11.7 

13.3 

14.7 

16.2 

17.7 

19.2 

20.6 

22.1 

22.4 

22.7 

22.9 

.000 

.002 

.015 

.017 

.018 

.016 

.021 

.029 

.034 

.033 

.022 

.012 

.004 

-.003 

-.008 

-.014 

-.022 

-.020 

-.022 

-.004 

.000 

256.0 

10.0 
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PROGRAM OUTPUT - SAMPLE CASE 


SPACE SHUTTLE BOOSTER REENTRY STANDARD ATMOSPHERE EASTERLY WINDS 


FLIGHT ALTITUDE = 50400, FEET 

MACH NUMBER = 1.200 

HEADING = 356.5 DEGREES 

FLIGHT PATH ANGLE = -12.75 DEGREES 

M-DOT = -.0197 /SEC 

PSI-DOT = -0.359 DEGREES/SEC 

GAMMA-DOT = 1.013 DEGREES/SEC 

AIRCRAFT LONGITUDE = 119.880 DEGREES 

AIRCRAFT LATITUDE = 27.600 DEGREES 

AIRCRAFT LENGTH = 256.0 FEET 


ATMOSPHERIC CONDITIONS 
POG = 2116.20 PSF 


ZT 

TO 

ZWE 

VOE 

ZWN 

VON 

FEET 

DEG F 

FEET 

FT/SEC 

FEET 

FT/SEC 

0.0 

59.0 

0.0 

5.0 

0.0 

0.0 

36200.0 

-69.7 

30000.0 

68.0 

202000.0 

0.0 

65799.9 

-69.7 

40000.0 

84.0 



105500.0 

-48.1 

45000.0 

79.0 



155500.0 

27.5 

55000.0 

36.0 



172000.0 

27.5 

65000.0 

19.0 



202000.0 

-4.8 

83000.0 

16.0 





110000.0 

34.0 
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5 • 30 


PHI 


47.00 DEGREES 


INITIAL WAVEFORM R/L = 


X/L 

DP/P 

0.000 

0.00000 

0.030 

0.00200 

0.060 

0.01500 

0.080 

0.01700 

0.100 

0.01800 

0.190 

0.01600 

0.230 

0.02100 

0.260 

0.02900 

0.280 

0.03400 

0.320 

0.03300 

0.490 

0.02200 

0.650 

0.01200 

0.790 

0.00400 

0.940 

-0.00300 

1.090 

-0.00800 

1.240 

-0.01400 

1.380 

-0.02200 

1.530 

-0.02900 

1.560 

-0.02200 

1.590 

-0.00400 

1 .610 

0.00000 


20 



WAVEFORM AT ALTITUDE 


30000.0 FEET 


LONGITUDE = 

119.922 

DEGREES 

LATITUDE = 

27.671 

DEGREES 

T, MSEC 


Pt PSF 

0.0 


0.000 

0.0 


1.477 

12.5 


1.412 

136.5 


0.770 

238.4 


0.257 

333.8 


-0.193 

411.7 


-0.513 

498.3 


-0.899 

551.1 


-1.164 

551.1 


0.000 

GROUND-RAY INTERSECTION 


LONGITUDE = 

120.051 

DEGREES 

LATITUDE = 

27.860 

DEGREES 

DISTANCE FROM GROUND 

TRACK = 9 


WAVEFORM AT THE GROUND REFLECTION FACTOR * 1.90 


T, MSEC P, PSF 


0.0 

0.0 

133.2 

258.6 

374.6 

467.1 
571.4 

608.2 
608.2 


0.000 
1.359 
0.789 
0.263 
-0. 197 
-0.525 
-0.919 
-1.074 
0.001 
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PROGRAM LISTING 


1 FORMAT ( 20AA ) 

2 FORMAT ( 10F7.0 ) 

3 FORMAT (315) 

A F0RMAT(////5X»2 AH SON IC CUTOFF ALTITUDE = ,F7.0,6H FEET, 

1 // 5X » 17HPH I FOR CUTOFF = ,F5.0,9H DEGREES) 

5 FORMAT ( 5X » 20AA, /// ) 

6 FORMAT (5X,18H FLIGHT ALTITUDE = ,F7.0,6H FEET , //5 X , 1 AHM ACH NUMBER 
1= ,F6.3,//5X,10HHEADING = ,F6.1,9H DEGR EES, //5X, 20 HF LIGHT PATH AN 
2GLE = , F6.2 »9H DEGREES, //5X,BHM— DOT = ,F6.A,6H /S EC , // 5 X , 1 OHPS I- 
3D0T = , F6 .3 , 1 3H DEGREE S/SEC , //5X , 12HGAMMA-D0T = ,F6.3,13H DEGREE 
AS/SEC, //5X»21HAIRCRAFT LONGITUDE = ,F8.3,9H DEGREES , //5X ,20HA I RCR 
5AFT LATITUDE = ,F7.3,9H DEGREES, //5X, 18HAIRCRAFT LENGTH = ,F5.1, 
66H FEET, ////5X,22HAT MO SPHERIC CONDITIONS) 

7 FORMATC/5X,6HPOG = ,F7.2,5H PSF) 

8 FORMAT (//9X,2HZT, 11X,2HT0» 11X, 3H ZWE , 1 1 X , 3HV0E, 11X,3HZWN, 11X,3HV0N/ 
18X,AHFEET ,9X »5HDEG F,9X , AHFEET ,9X , 6HF T/SEC ,9X , AHFE ET , 9X , 6HFT/ SEC/ ) 

9 FORMAT (F1A.1,F11.1»F15.1»F12.1»F16. 1,F12. 1) 

10 F0RMAT(25X»F15.1»F12.1»F16.1»F12. 1 ) 

11 FORMAT (25X,F15. 1,F12. 1) 

12 F0RMAT(52X,F16.1,F12.1) 

13 FORMAT (F1A.1,F11.1,27X,F16.1,F12.1) 

1A F0RMAT(F1A,1 , FI 1 . 1 ) 

15 FORMAT (F1A.1,F11.1,F15.1,F12.1) 

16 FORMAT (//////5X,16HINITIAL WAVE FORM , 1 3X , 6HR/L = , F5. 2, 1 3X,6HPH I = 

1 , F7. 2 , 8H DEGREES, /// 1 6X , 3HX /L , 17X, AHDP/P, //< F20. 3, F2 1.5) ) 

19 FORMAT ( 1 HI ) 

20 FORMAT ( ///5X » 23HWAVE FORM AT ALTITUDE = ,F8.1,6H FE ET, // 9X , 1 2HL0N 
1GITUDE = ,F8.3,9H DEGREES, /9X , 11HLATI TUDE = ,F7.3,9H DEGREES,//) 

21 FORMAT I 13X,7HT , MSEC , 15X , 6HP , PSF , // ( F20 . 1 , F20 . 3 ) ) 

22 FORMAT (13X,7HX, FE ET , 14X , 6HP , PSF , // ( F20 . 2 , F20. 3 ) ) 

23 FORMAT ( 13X»7HT , MSEC , 15X , AHDP/P, //(F20.1,F20.6)) 

2 A FORMAT (13X»7HX» FEET , 1AX , AHDP /P , / / ( F20. 2 , F 20 . 6 ) ) 

27 FORMAT (8X ,3AHR AY TUBE AREA GOES TO ZERO AT Z = ,F8.1,6H FEET) 

28 FORMAT (/5X,28HC0NSTANT PRESSURE ATMOSPHERE) 

29 FORMAT ( ///5X , 23HGR0UND— RAY I NT ERS ECT I ON , //9 X , 1 2HL0NG I TUD E = ,F8.3 
1 , 9H DEGREES ,/9X ,11HLAT ITUDE = ,F7.3,9H DEGR EES , /9X , 29HD I S T ANC E F 
2 ROM GROUND TRACK = ,F6.2,7H M I L ES, //// 5X , 22HWA VE FORM AT THE GROUN 
3D,10X»20HREFLECTI0N FACTOR = ,FA.2,//) 

INTEGER OCODE 

REAL MACH,NI3,A),K1,M( 100),LAMDA< 100 ) , ML, MAG, MOOT , LATP , LONGP 
DIMENSION TITLEI20) ,R( 3, A) »ZT( 100) ,T0( 100) , ZWE ( 100) ,VOE( 100) , 

1 ZWN( 100) ,VON( 100 ) , ALT (50) ,X( 100 ) ,DPP( 100) , FI ( 100) , F2( 100 ) , 

2 P< 100) ,DP ( 100) ,XX( 100 ) ,PP ( 100) ,V0( 3) , A01( A) , V01 ( 3, A) 

DOUBLE PRECISION MM ACH , FFPA , HH EAD, MU 

WRITE ( 6 , 19 ) 

AO READ (5,1) TITLE 

Rf^ADI 5,2) MACH, FLT ALT 

READ (5,2) MDOT , PSIDOT, GAMDOT, FPA 

READ (5,2) LONGP, LATP, HEAD 

READ(5»2) POG 

READ (5,3) NTEMP, NWINDE, NWINDN 

RE AD ( 5 ,2 ) (ZT(I), 1=1, NTEMP) 

READ! 5,2 ) <TO(I), 1 = 1, NTEMP) 

RE AD ( 5 , 2 ) ( ZWE ( I ) , 1=1, NWINDE) 

READ (5,2) (VOE(I), 1=1, NWINDE) 

RE AD ( 5 , 2 ) (ZWN(T), 1=1, NWINDN) 

READ (5,2) (VON(I), 1=1, NWINDN) 

READ ( 5 ,3 ) NALT 

READ (5,2) ( ALT ( I ) , 1 = 1, NALT) 
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45 


50 

60 

70 

80 


85 

90 

100 


110 

120 

130 

140 

150 

160 

170 

180 


READ ( 5 , 3 ) 
R EAD ( 5 * 2 ) 
RE AD ( 5 , 2 ) 
R EAD ( 5 ,3 ) 
READ (5*2) 
READ ( 5 , 2 ) 


OCODE 

REEL 

ROVERL, PHI 
NX 

( X ( I ) ♦ I = 1 , NX ) 
(DPP(I)t 1 = 1 , NX ) 


RE AD ( 5 ,2 ) ALt ML 

DELTAT=. 1 
PH I S AV=PH I 
R0=R0VERL#A L 

G= 32 • 2 6-0# 17*C0S(LATP/57.296)**2 
I F ( POG .GT. 0.) GO TO 45 


G = 0. 
POG = 1 • 


ICUT=0 


DPH I = - 10 • 

IF (PHI .LT. On DPH 1 = 10* 
XOFF = X(l) 

DO 50 1=1, NX 

X( I ) = ( X C I )-XOFF)/ML 


DO 60 1=1 ,NTEMP 

ZT ( I ) = 1 000 • 0*ZT ( I ) 

DO 70 I =1 ,NW I NDE 
ZWE ( I ) =1000 • 0*ZWE ( I ) 

DO 80 1 = 1 , N W I NON 

ZWN< I ) =1000 .0*ZWN ( I ) 

WR I TE ( 6 , 5 ) TITLE 

WRITE (6,6) FLTALT,N!ACH,HEAD,FPA,MDOT,PSIDOT,GAMDOT,LONGP,LATP,AL 
I F ( G .GT. 0. ) GO TO 85 
WRITE (6,28) 

GO TO 90 

WRITE (6t7) POG 
WRITE(6,8 ) 

1 = 1 

I F ( I .GT. NTEMP ) GO TO 110 
I F ( I .GT. NWINDE) GO TO 140 

I F ( I .GT. NWINDN) GO TO 160 

WRITE (6,9 ) ZT(I ),T0( I ) ,ZWE (I ) , VOE ( I ) , ZWN( I ) , VON( I ) 

GO TO 170 

I F ( I .GT. NWINDE) GO TO 130 

I F ( I .GT. NWINDN) GO TO 120 

WRITE (6, 10) ZWE (I ),VOE( I ) , ZWN ( I ) , VON ( I ) 

GO TO 170 

WRITE (6,11 ) ZWE ( I ), V0E( I ) 

GO TO 170 

I F ( I .GT. NWINDN) GO TO 180 
WRITE (6,12) ZWN ( I ) ,VON( I ) 

GO TO 170 

I F ( I .GT. NWINDN) GO TO 150 

WRITE (6,13) ZT ( I ) ,T0( I ) , ZWN (I ) , VON ( I ) 

GO TO 170 

WRITE (6,14) ZT ( I ),T0(I) 

GO TO 170 

WRITE (6, 15) ZT(I ),TO(I),ZWE( I ) , VOE ( I ) 

1 = 1 + 1 
GO TO 100 

I F ( I .LT. 10) WRITE (6,19) 

WRITE (6,16) ROVERL ,PH I , ( X ( I ) , DPP ( I ) , I = 1 , NX ) 

WR ITE ( 6 , 19 ) 
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200 I F ( I CUT .EG. 0) GO TO 230 

I F ( PH I S A V . LT. 0. ) GO TO 210 
I F ( PHI .GT. 0. ) GO TO 230 
DPH I =1 • 

GO TO 220 

210 IF ( PHI .LT. 0.) GO TO 230 
DPH I =— 1 • 

220 PH 1 = 0. 

230 MMACH=MACH 

FFPA=FPA/57.296 
HHEAD=HEAD/57.296 
PHI=PHI/57.296 
R ( 1 , 1 ) =0 • 

R ( 2 1 1 ) =0. 

R ( 3 , 1 ) =FLTALT 
I ZER0 = 0 
I FORK= 1 
11=1 
I R AY=2 
H= FLTALT 

240 DO 250 J = 2 , NT EMP 

IFCZT(J) .LT. H) GO TO 250 

A0=49.1*SQRT(T0( J-1)+(H-ZT( J-l) )/(ZT( J)-Z T( J- 1 ) > * ( T0( J ) -TO ( J- 1 > ) 
1+459.67) 

GO TO 260 
250 CONTINUE 
260 DO 270 J=2,NWINDE 

IF ( Z WE ( J ) .LT. H) GO TO 270 

VO ( 1)=V0E( J-1)+(H-ZWE( J-l) ) /(ZWE( J)-ZWE( J-l) ) * < VOE ( J ) -VOE ( J- 1 ) ) 
GO TO 280 
270 CONTINUE 
280 DO 290 J = 2 ,NW INDN 

I F ( Z WN ( J ) .LT. H) GO TO 290 

VO ( 2 ) =V0N( J-l ) + ( H-ZWN ( J-l ) ) /( ZWN ( J )-ZWN ( J-l ) )*( VON( J ) -VON ( J- 1 ) ) 
GO TO 300 
290 CONTINUE 

300 I F ( I RAY .GT. 1 .OR. ICUT .EQ. 1) GO TO 340 
P0= POG 

DO 330 J = 2 * NT EMP 
Z2=ZT ( J ) 

IF ( H .LE. Z 2 ) Z 2 = H 

I F { TO ( J-l ) .EO. TO ( J ) ) GO TO 310 
TOZ=(TO( J)-TO( J-l ) ) /( ZT( J)-ZT( J-l) ) 

PO=PO*( l.+TOZ/(TO( J-l) +459.67 )*(Z2-ZT( J-l ) ) ) ** ( - 1 . / 53 . 3 5 3/TO Z ) 
GO TO 320 

310 PO=PO*EXP { ( ZT ( J-l ) -Z2 ) / 53 • 35 3/ ( TO ( J-l ) +459 .67 ) ) 

320 I F ( Z 2 .EO. H) GO TO 340 
330 CONTINUE 

340 GO TO (350,460,580), I FORK 
350 VOXH= VO ( 1 ) 

VO YH= VO ( 2 ) 

VX=MACH*A0*C0S(FPA/57.296)*SIN(HEAD/57.296 ) 
VY=MACH*A0*C0S(FPA/57.296)*C0$(HEAD/57.296) 

IR A Y=0 

360 IRAY= I RAY+ 1 

GO TO (410,390,380,370,430), IRAY 
370 R(1,4)=R(1,3) 

R ( 2 1 4) =R ( 2 , 3 ) 

R(3,4)=R(3,3) 

GO TO 400 
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380 


DT=0.17453*RO/AO/SQRT(MACH**2-1.0) 

MMACH=MMACH+MDOT*DT 
FFPA=FFPA+GAMD0T/57.296*DT 
HHEAD=HHEAD+PSIDOT /57. 2 96*01 
PHI=PHI-0. 017453 
R< 1 ,3)=R ( 1 ,1 )+(VX+V0( 1 ) )*DT 
R(2,3)=R(2,1 >+(VY+V0( 2) )*DT 
R (3,3 )=R (3 ,1 )+MACH*A0*SIN(FPA/57.296)*DT 
GO TO 410 
390 R(1,2)=R(1»1) 

RC2»2)=R(2,1) 

R(3,2)=R<3,1) 

400 PHI=PHI+0. 017453 
GO TO 420 

410 MU=DARSIN( 1.0/MMACH) 

420 N( 1, IRAY)=DCOS< FF P A ) *DSI N ( HHE AD ) *DS I N( MU ) -DCOS ( HHEAD )*DCOS ( MU ) 

1*S I N ( PH I )+DSIN(FFPA)*DSIN( HHEAD)* DCOS( MU )*COS( PHI ) 

N ( 2 * I RAY ) =DCOS ( FFPA)*DCOS(HHEAD)*DSIN(MU)+DSIN<HHEAD)*DCOS(MU) 

1 *S I N ( PH I )+DSIN(FFPA)*DCOS(HHEAD)*DCOS(MU)*COS< PHI ) 

N(3, IRAY)=DSIN< FFPA)*DSIN(MU)-DCOS( FFPA)*DCOS(MU)*COS(PHI ) 

I F ( I RAY .EC. 4) PHI=(PHI-.017453)*57.296 

GO TO 360 

430 DO 440 I R AY= 1.4 

I F ( I RAY .EQ. 1) DELT=R0/A0/0C0S( MU) 

I F ( I RAY .EC. 3) DELT =DELT-DT 
R ( 1 » I RAY ) =R ( 1 , IRAY) + ( AO*N( 1, IRAY) + VO( 1) >*DELT 
R( 2 ,IRAY )=R (2 ,IRAY )+( A0*N( 2,IRAY)+V0i 2) )*DEi_T 
440 R(3, IRAY)=R<3, IRAY )+AO*N ( 3, IRAY )*DFLT 

Al=( (R (2,4)-R(2, 1 ) )*(R(3,3)-RC3,2) )-(R(3,4)-R( 3, 1) )*(R( 2,3) 
l-R(2,2)))*N(l,l)+< (R(3,4)-R(3, 1>)*(R(1,3)-R(1,2))-(R(1,4)-R<1,1))* 
2 ( R ( 3 ,3 )— R ( 3 » 2 ) ) ) *N ( 2 * 1 ) + ( ( R ( 1, 4) -R C 1 » 1 ) )* ( R ( 2 , 3 ) — R ( 2, 2 ) ) 
3-<R(2,4)-R(2,l ) >* (R( 1,3)-R ( 1,2 ) ) )*NI 3, 1 ) 

I F0RK=2 
I R A Y = 0 

450 IRAY=IRAY+1 
H=R ( 3 » IRAY ) 

GO TO 240 
460 AO 1 ( IRAY ) = AO 

V01( 1, IRAY)=VO( 1 ) 

V01(2»IRAY)=V0(2) 

I F ( IRAY .LT. 4) GO TO 450 

I F ( I CUT .EQ. 1) GO TO 530 

P01=P0 

CN1=A01(1)+V01(1,1)*N<1,1)+V01(2,1)*N(2,1) 

DO 470 1 = 1, NX 

P( I ) =DPP ( I )*P0 
470 X( I )=X(I )/MACH*AL/CNl 
K= 1 
J=2 

NN=NX-1 
MCI) =0. 

DO 520 1=1, NX 

I F ( K .EQ. 2) GO TO 510 
I F ( I .EQ. NX ) GO TO 500 
I F ( X ( I > .EQ. X ( I + 1 ) ) GO TO 480 
DP ( J ) =0 . 

LAMDA ( J )=X ( 1+1 ) — X ( I ) 

M(J) = (P(I+1)-P( I ) ) /L AMD A( J ) 

J=J+1 
GO TO 520 
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480 K = 2 

DP( J)=P(I+1 )-P( I ) 

I F C 1+2 ♦GT. NX) GO TO 490 
LAMDA (J)=X(I+2)-X(I+l) 

M( J) = (P ( 1+2 )-P( 1 + 1 ) ) /LAMDA ( J ) 

NN=NN- 1 
J = J + 1 
GO TO 520 
490 M { J ) =0 • 

NN=NN- 1 
GO TO 520 
500 DP ( J ) =0 • 

M ( J ) =0* 

510 K= 1 
520 CONTINUE 
530 KR I TE=0 

I F ( DELTAT .EQ. • 2 ) GO TO 535 

DELTAT^. 0002* SORT ( R ( 1,1 )**2+R ( 2, 1 ) **2 + ( FLTALT-R ( 3, 1 ) )**2) 

I F ( DELTAT .GT. .2) DELTAT= . 2 
53 5 I F ( R ( 3 , 1 ) + AO 1 ( 1 ) *N ( 3 , 1 > *DELT AT-ALT ( II ) ) 5 40, 540,550 
540 DELT=(R (3,1 >-ALT( II ) ) /AO 1 ( 1 ) / ( -N ( 3, 1 ) ) 

KRI TE= 1 
GO TO 560 
550 DELT=DELTAT 
560 I RAY =0 

AOS = AO 1 ( 1 ) 

570 IRAY= I R AY+1 

I F ( IRAY .EO. 5) GO TO 610 

R(1 , I RAY ) =R ( 1,IRAY ) + (A01( IRAY)*N( 1,IRAY)+V01( 1,IRAY) ) *DELT 
R ( 2 , IRAY)=R(2, IRAY)+( A01 ( I R AY ) *N ( 2 , IRAY)+V01( 2, IRAY) )*DELT 
R(3,IRAY )=R(3,IRAY J+AOl ( I R AY ) *N ( 3 , 1 R AY )*DELT 
H=R ( 3 , I R AY ) 

I E0RK = 3 
GO TO 240 

580 DZ=-A01 ( IRAY )*N(3, IRAY ) *DELT 

AOZ= ( A01 ( I R AY ) — AO ) /DZ 
VOXZ=( V01 ( 1 , IRAY )-V0( 1 ) ) /DZ 
VOYZ=(V01(2,IRAY)-V0( 2)) /DZ 
A01< IRAY )=A0 
VO 1(1, IRAY) =V 0(1) 

V01(2,IRAY)=V0(2 ) 

CC=- ( N ( 1 , IRAY ) *VOXZ+N ( 2, IRAY ) *V0YZ+A0Z ) 

DN1=-N(1 ,IRAY)*N(3,IRAY)*CC*DELT 
DN2=-N(2, IRAY)*N(3,IRAY)*CC*0ELT 
DN3=(N( 1 tIRAY )**2+N(2, IRAY ) **2 ) *CC*DELT 

MAG=SQRT ( (N( 1, IRAY )+DNl)**2+(N( 2,IRAY)+DN2)**2+( N( 3, IRA Y ) +DN3 ) **2 ) 
N(1,IRAY )=(N( 1,IRAY )+DNl )/MAG 
N ( 2 , IRAY) = (N(2, IRAY )+DN2) /MAG 
N(3,IRAY )= (N(3, IRAY )+DN3) /MAG 
IF (N ( 3 , IRAY ) .LT. -.01745) GO TO 570 
I F ( ICUT .EO. 0) GO TO 600 
IF ( ABS ( DPH I ) .GT. 2.) GO TO 590 
WR I T E ( 6 , 4 > HCUTtPHI 

GO TO 1090 
590 PHI=PHI+DPHI 
GO TO 200 
600 HCUT=R (3,1) 

ICUT= 1 
I PH I=PH I 
PH I - I PH I 
GO TO 590 
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610 IFUCUT .EQ. 1 .OR. I ZERO .EO. 1) GO TO 921 
CN =A01 ( 1 ) +V01 ( 1*1 )*N( 1*1)+V01(2, 1 )*N ( 2, 1) 

A = ( (R(2»4)-R(2,l ) )*(R(3,3)-R<3,2) )-(R(3,4)-R(3,l) )*(R(2,3) 
1-R(2,2)) )*N(l,l)+( (R(3,4)-R(3,1))*(R( 1»3)-R( 1»2) )-<R< 1»4)-R( 1,1) )* 
2(R<3,3 )-R(3,2) ) )*N(2, l) + ( (R(l,4)-R(l, 1) )* (R ( 2, 3 )-R( 2, 2 ) ) 
3-(R(2,4)-R(2,l))*(R(l»3)— R(1»2)))*N(3»1) 

I E ( A .GT. 0.) GO TO 620 
HZERO=R (3,1) 

I ZERO* 1 
GO TO 921 

620 Cl* 1.71 4* (AO 1(1) +AOS )/(P01+P0) /( CN+CN1 ) 

C2= ( ( A01 ( 1 )-A0S ) /( A01 ( 1) +AOS) +2.8*DZ*G/( A01 ( 1 ) + AOS ) **2-2.* ( CN-CN1 ) 
1/ ( CN+CN1 )— (A— Al)/( A+Al ) J/OELT 
P01=P0 
CN1=CN 
A 1=A 

630 DELTS=DELT 

TP=(EXP(C2*DELT)-1.0)/C2 
N1=NN+1 
N2=NN+2 
K = N2 
640 Jl=l 

650 00 660 J=1,N2 

I F ( M ( J ) .EQ. M(K) .AND. J1 .EO. 2 .AND. J .NE. K) M ( J ) = .9999*M ( J ) 
Fl( J)=1.0-C1*M( J)*TP 

IF(F1(J) .LE. 0. .AND. J .NE. K) GO TO 670 
660 CONTINUE 
GO TO 680 
670 K= J 
J 1 = 2 

TP=1.0/C1/M( J ) 

GO TO 650 
680 TK=C1*TP 

IF ( J 1 .EQ. 1) GO TO 690 

I F ( DP ( K ) .EQ. 0. .AND. 0P(K+1) .EQ. 0.) GO TO 690 
IF ( DP ( K ) .EQ. 0.) GO TO 760 
I F ( DP ( K+l ) .EQ. 0.) GO TO 780 
GO TO 800 
690 KF= 1 
J=1 

695 J*J+1 

IF ( J .GT. Nl) GO TO 840 
I F( J .EQ. K) GO TO 695 

700 IF(ABS(TK*(M( J)-M( J-l) l/Fll J) ) .GT. .001) GO TO 710 

IF ( ABS (TK*(M( J )-M (J+1))/F1(J) ) .GT. .001) GO TO 705 

F2( J)=LAMDA( J )-(DP( J)+DP( J+l) )*TK/F1( J) /2. 

GO TO ( 730*820 ) * KF 

705 F2( J)=LAMDA( J)-OP( J)*TK/F1( J)/2.-0P( J+l)/( M( J)-M( J+l) ) * ( SORT ( FI ( J+ 

1 1 ) /FI ( J ) ) — 1. ) 

GO TO (730,820), KF 

710 IF ( ABS ( TK* ( M ( J )-M ( J+l ) ) /FI ( J ) ) .GT. .001) GO TO 720 

F2 ( J ) =LAMDA ( J )-DP ( J ) /( M ( J )-M ( J-l ) )* ( SQRT( F 1 ( J-l ) /FI ( J ) ) -1 . ) -DP( J+l 
1)*TK/F1( J)/2. 

GO TO (730,820), KF 

720 F2(J)=LAMDA(J) -DP ( J ) / ( M( J )— M( J-l ) )* ( SORT ( FI ( J-l ) /FI ( J ) )-l . ) — OP ( J+l 
1)/(M( J)— M( J+l ) )*(SQRT(F1( J+ll/Fll J))-l. ) 

GO TO (730,820), KF 
730 IF(F2(J)> 740, 735, 695 
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735 I F ( J 1 .EQ. 2) L AMDA ( J ) = 1 . 000 1*LAMDA ( J ) 

I F ( J 1 .EQ. 2) GO TO 700 
K = J 

GO TO 695 
740 K= J 

750 I F ( DP ( K ) .GT. 0.) GO TO 770 

IF(ABS(TK*(M(K)— M(K+1))/F1(K) ) .GT. .001) GO TO 760 
TP=2.*LAMDA(K > /Cl / < DP ( K+l ) +2 . *M < K )*LAMDA( K) ) 

GO TO 640 

760 TP=( 1.— (1. +LA MDA (K)/DP(K+1)*(M(K) — M ( K+ 1 ) ) ) ** 2 ) / ( C1=M M ( K+ 1 ) — MC K ) * 

1 ( l.+LAMDA(K ) /DP (K+l )* (M ( K) -Ml K+l ) ) ) **2 ) ) 

GO TO 640 

770 IF C DP (K + l ) .GT. 0.) GO TO 790 

IF(ABS(TK*(M(K)-M(K-1) ) /F 1 ( K » ) .GT. .001) GO TO 780 
TP=2.*LAMDA(K)/C1/(DP(K)+2.*M( K ) *L AMDA ( K) ) 

GO TO 640 

780 TP=(1.-<1.+LAMDA(K)/DP(K)*(M(K)-M<K-1) ) ) **2 ) / ( C 1* ( M ( K- 1 ) -M ( K ) * ( 1 . 

1+LAMDA(K)/DP(K)*(M(K)— M(K-l) ) )**2) ) 

GO TO 640 

790 IF(ABS(TK=MM(K)-M(K-1) )/Fl(K) ) .GT. .001) GO TO 800 

IF(ABS(TK*(M(K)-M(K+1))/F1(K> ) .GT. .001) GO TO 800 

TP=2.*LAMDA(K ) /Cl / ( DP ( K ) +DP ( K+l ) +2. *M I K )*LAMDA I K) ) 

GO TO 640 
800 TP1=0. 

TP2=TP 
KF = 2 
IC = 0 

810 IC=IC+1 

I F ( IC .GT. 20) GO TO 640 
TP= ( TP1+TP2 )/2. 

TK=C1*TP 

F1(J)=1.0-C1*M(J)*TP 
F 1 ( J- 1 ) = 1 . 0-C 1#M ( J— 1 ) *T P 
F 1 ( J + l ) =1 ,0-C 1 *M ( J+ 1 ) *TP 
GO TO 700 

820 I F ( F2 ( K ) .LE. 0.) GO TO 830 
TP1=TP 
GO TO 810 
830 TP2=TP 

GO TO 810 

840 DELT=AL0G(TP*C2+1.0)/C2 

E2=EXP (C2*DELT) 

I F ( K .EQ. N2 ) GO TO 850 
DPSUM=DP(K)+DP(K+1) 

I F ( K .EQ. 2) GO TO 870 
850 K 1 1=K-1 

DO 860 J = 2 »K 1 1 
M<J)=M(J)/F1(J)*E2 

DP ( J ) =DP ( J )*E 2/ SORT (F1(J)*F1(J— 1) ) 

860 LA MDA (J)=F1(J)*F2(J) 

I F ( K .EQ. N2 ) GO TO 920 
870 K22=K+1 

I F ( DPSUM .EQ. 0.) GO TO 880 

DP ( K ) =E2* (DP(K )/SQRT (F1(K)*F1(K-1) ) +DP ( K+ 1 ) /SQR T( F1(K)*F1(K+1))) 
GO TO 890 

880 DP(K)=M(K)*tAMDA(K)*E2 
890 M ( K )=M ( K+l ) /F 1 ( K+l )*E2 

I F ( K .EQ. Nl) GO TO 910 

DO 900 J=K22»N1 

M( J)=M( J+1)/F1( J+1)*E2 

DPI J)=DP( J + l )*E2/SQRT( F 1 ( J+l )*F1( J) ) 
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900 LAMDA(J-1)=F1(J)*F2( J) 

910 NN=NN-1 

DELT=DELTS-DELT 
GO TO 630 

920 DP(N2)=DP(N2)*E2/SQRT(F1(N2)*F1(N1) ) 

921 I F ( KR I TE .EO. 0) GO TO 530 
I F ( ICUT .EO. 0) GO TO 922 
11 = 11+1 

I F ( II .LE. NALT) GO TO 530 
OPH I = 1 . 

I F ( PHI SA V .LT. 0.) DPHI =-l • 

GO TO 590 

922 I F ( I ZERO .EQ. 0) GO TO 926 
11 = 11+1 

I F ( II .LE. NALT) GO TO 530 
I F ( ALT ( NALT ) .NE. 0.) GO TO 926 
11 = 11-1 
GO TO 990 
9?6 NALT = NALT + 1 
ALT ( NALT ) =0. 

GO TO 530 
926 XX(1)=0. 

PP( 1 )=0. 

J = 2 

DO 960 L=2,N1 

I F { DP ( L ) .GT. 0.) GO TO 930 
XX ( J)=XX( J-1)+LAM0A(L) 

PP< J)=PP< J-1)+M(L)*LAMDA(L) 

J=J+1 
GO TO 960 
930 XX( J)=XX(J-1) 

PP< J)=PP( J- 1 » +DP ( L ) 

XX ( J+1)=XX (J)+LAMDA( L) 

PP(J+1)=PP(J)+M(L)*LAMDA(L) 

J= J + 2 

960 CONTINUE 

I F ( DP ( N2 ) .EQ. 0.) GO TO 950 
XX(J)=XX(J-1) 

PP(J)=PP( J-1)+DP(N2) 

J=J + 1 
950 J=J-1 

I F ( ALT (II) .NE. 0.) GO TO 970 
DO 960 L = 1,J 
960 PP(L)=PP(L)*REFL 
970 I F ( POG .NE. 1.) GO TO 990 
DO 980 L=1 , J 

980 PP(L)=PP(L)/PO 

990 RLAT=R (2*1 )/( 20890070. + ALT ( 1 1 ) )*57.296+LATP 

RLONG=-R (l»l)/(20890070.+ALT( II ) >*57. 296/COS ( RLAT/57. 296 ) +LONGP 
IF ( ALT (II) .NE. 0.) GO TO 1000 

DGT= (R (2,1 )*SIN(HEAD/57.296)-R( 1 , 1 ) *COS ( HEAD/57. 296 ) J/5280. 
WRITE(6,29) RLONG,RLAT,DGT,REFL 
I F ( I ZERO ) 1010, 1010, 1080 

1000 WRI TE( 6 , 20 ) ALT ( 1 1 ) ,RLONG, RLAT 
1010 I F ( OCODE .EQ. 2) GO TO 1060 
DO 1020 L = 1 , J 
1020 XX (L ) =XX ( L ) *1000. 

I F ( POG .EQ. 1.) GO TO 1030 
WRITE (6,21 ) (XX(L),PP(L)»L=1»J> 

GO TO 1070 


1030 WRITE (6*23) < XX ( L ) , PP ( L) , L= 1 , J ) 

GO TO 1070 

1040 C= ( ( VX+VOXH ) *N ( 1 , 1 ) + ( VY+VOYH ) *N ( 2, 1 ) ) /SORT! ( VX+VOXH) **2+ 
1(VY+V0YH)**2) 

DO 1050 L = 1 » J 
1050 XX(L)=XX<l)*CN/C 

I F ( POG .EQ. 1. ) GO TO 1060 
WRITE(6»22) (XX(L)tPP(L) *L=1»J) 

GO TO 1070 

1060 WRI TE ( 6 * 24 ) ( XX ( L ) , PP ( L ) , L= 1 * J ) 

1070 11=11+1 

KRITE=0 

IF(II-NALT) 530, 530, 1090 
1080 WRITE(6»27) HZERO 

1090 WR ITE ( 6 , 19 ) 

GO TO 40 
END 
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